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1. Introduction 

The solar magnetic field is an important quantity which couples the so- 
lar interior, the photosphere and the atmosphere. The quasi stationary 
coronal magnetic field configuration is an interesting and challenging 
topic on its own right. But even to understand basic processes like coro- 
nal mass ejections and flares it is important to understand the quiescent 
magnetic configuration out of which these dynamic phenomena arise. 
Unfortunately the coronal magnetic field cannot be measured directly, 
but it has to be reconstructed from photospheric measurements. A mag- 
netic field reconstruction of the solar corona has to be consistent with 
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the observed spatial variation of the coronal plasma (density, pressure, 
temperature) often elongated along the magnetic field. 

Here we are mainly interested in long living structures which are 
time independent in first order. We also concentrate on closed magnetic 
configuration where a stationary plasma flow (solarwind) does not sig- 
nificantly contribute to the force balance. Such configurations are static 
equilibria and have to obey the magnetohydrostatic equations (MHS). 

As the magnetic field B and the density distribution ./V are physically 
closely related their model reconstruction should also be linked as much 
as possible. In this paper we attempt to show how this can be achieved. 
We propose variational principles which if they can be solved should 
give a consistent model for an isothermal corona. For the magnetic field 
reconstruction this leads to a generalization of a nonlinear force-free 
approach by (Wheatland et al., 2000), for the density reconstruction 
we obtain a tomography problem with an improved regularization term. 

The ground based or space-born magnetograph observation provide 
either the line-of-sight magnetic field (Bi os , e.g. MDI on SOHO), which 
is sufficient for potential and linear force-free fields, or all three compo- 
nents of the photospheric magnetic field (e.g. IVM in Hawaii, expected 
also from SolarB). The latter information is sufficient to determine 
nonlinear force-free fields completely. As a force-free approximation is 
justified only in the limit of a vanishing plasma (3, we take into account 
forces (pressure gradient and gravity) for configurations with a finite 
plasma (5 even though we shall consider (3 small. 

Popular simplifications for the reconstruction of coronal magnetic 
fields are: 

— Potential fields (j = 0) (e.g. Schmidt, 1964; Semel, 1967; Schatten, 
Wilcox, and Ness, 1969; Sakurai, 1982; Rudenko, 2001a) 

— Linear force-free fields (e.g. Nakagawa and Raadu, 1972; Chiu and 
Hilton, 1977; Seehafer, 1978; Semel, 1988; Gary, 1989; Lothian and 
Browning, 1995) 

— Linear non force-free fields (e.g. Zhao and Hoeksema, 1993; Zhao 
and Hoeksema, 1994; Petrie and Neukirch, 2000; Zhao, Hoeksema, 
and Scherrer, 2000; Rudenko, 2001b) 

— Nonlinear force-free fields (e.g. Sakurai, 1981; Wu, Chang, and 
Hagyard, 1985; Roumeliotis, 1996; Amari et al., 1997; McClymont, 
Jiao, and Mikic, 1997; Wheatland et al., 2000; Yan and Sakurai, 
2000). 

Within this work we do not use any of these assumptions but consider 
the general case of nonlinear non-force-free equilibria. The mathemat- 
ical problem of calculating nonlinear non-force-free fields is closely 
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related to the problem of calculating nonlinear force-free fields which 
coincides with the above in the limit of (3 — ► 0. Under ideal con- 
ditions the information contained in a (perfect) vectormagnetogram 
together with the force-free condition would be sufficient to calculate 
the coronal magnetic field. Within this work we show that the informa- 
tion contained in a vectormagnetogram together with a tomographic 
reconstructed coronal density distribution and the assumption of mag- 
netohydrostatic force balance is as well sufficient to calculate the finite (5 
coronal magnetic field. Unfortunately current vectormagnetograms and 
tomographic reconstruction are far from being perfect, which affects the 
quality of reconstruction. Within this work we use well known MHS- 
equilibria to test our newly developed reconstruction program. The 
use of analytic equilibria as artificial data allows us to extract ideal 
vectormagnetograms as well as ideal coronal density distributions. 

As for the density observations, ground based coronagraphs (e.g., 
the Mark III coronagraph on Hawaii, LASCO coronagraph on SOHO 
and the future STEREO mission) provide the line of sight integrated 
density structure of the solar corona from different relative viewpoints 
as the Sun rotates. These measurements have been used for a 3D- 
reconstruction of the coronal plasma distribution with help of tomo- 
graphic methods (Davila, 1994; Zidowitz, 1999; Frazin, 2000; Frazin, 
2002). The major problems here are: 

— the assumption of stationarity of coronal structures as the Sun 
rotates, 

— the lack of data due to the occulted center of the image, 

— the nonideal viewing geometry caused by a slight tilt of the Sun's 
axis with respect to the ecliptic. 

These shortcomings generally enhance the intrinsic ill-posedness of the 
tomography problem. The general approach to stabilize the reconstruc- 
tion is to smooth the solution by regularization. The prize to pay is a 
reduced spatial resolution of the model depending on the quality of the 
data and inconsistencies and ill-conditioning due to the above effects. 
So far only very general, isotropic regularization operators have been 
applied to coronal density reconstruction problems. Our approach to 
the density reconstruction in connection with the reconstruction of the 
coronal magnetic field leads to a new regularization operator which, 
as we demonstrate by test calculations, could yield a better spatial 
resolution than conventional reconstructions. 

The paper is outlined as follows. In section 2 we describe the basic 
equations and the newly developed algorithm of the magnetic field re- 
construction program in the case where the plasma density distribution 
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N is given. Section 3 contains several test-runs where we apply our 
code for the reconstruction of analytic MHS-equilibria. In section 4 
we propose an algorithm for an improved reconstruction of N if some 
information of B is given. This approach is tested and compared with 
conventional methods by with the help of a two dimensional analytic 
coronal density distribution. In the final section 5 we discuss how both 
methods could be used together to derive a consistent model of the 
Sun's corona. In appendix A, B and C we provide the algebra which 
has been omitted in the text. 



We describe the coronal plasma with help of the magnetohydro static 
(MHS) equations. The MHS equations are 



where B is the magnetic field, j the electric current density, P the 
plasma pressure, p the plasma density, (Iq the vacuum permeability 
and ^ the solar gravity potential. We define the functional 



The domain V is a volume which on one side is bounded by the 
sun's photosphere. Obviously, L is bound from below by 0. This bound 
is attained if the magnetic field satisfies the MHS equations. Here we as- 
sume that the plasma pressure and the density are given. It is assumed 
that the corresponding information will be provided by tomographic 
reconstruction of the solar corona. We vary functional L with respect 
to an iteration parameter t and get (see Appendix A for the derivation) 



2. Basic equations 



j x B - VP - pW = 0, 
V x B = /ioj, 
V B = 0, 



(1) 
(2) 
(3) 






where 



F = Vx(fi a xB)-n a x(VxB) 

+V(ft b • B) - n b (v • B) + (fi£ + ng) B, 



(6) 



G = n x (O a x B) - n(ft b • B) 



(7) 
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ft a = B- 2 [(V x B) x B - W (VP + pV*)] , 

« b = B- 2 [(V B) B] (8) 

The surface integral in (5) vanishes if B is prescribed on the boundary. 
We iterate the magnetic field inside the computational box with 

which ensures that L is monotonously decreasing. For the bottom the 
boundary values are given by the photospheric vector magnetograph 
observations. On other boundaries we may either assume B or include 
the boundary values in the variation. Actually the handling of the not 
observed lateral and top boundaries of the computational box is similar 
here as in the nonlinear force-free case (Wiegelmann and Neukirch, 
2003). On the boundary of the computational box the magnetic field is 
iterated with 

<9B 

-7^- = where B observed, (10) 
<9B 

— = n G else. (11) 

We propose to use this iteration process to solve for the minimum 
of L. If a solution of the MHS-equations for the prescribed boundary 
condition exist, the global minimum of L corresponds to this solution 
and attains L = 0. Please note that the iteration procedure ensures 
to find this global minimum if the solution space is convex. For a 
non convex solution space it is possible that the iteration will lead 
to a local minimum. For complicated magnetic field configurations it 
is difficult to decide in advance whether the solution space is convex 
or not. For a non convex solution space it is still possible to find the 
global minimum by iteration if the start configuration is sufficient close 
(within a convex area) to this minimum. The method generalizes an 
approach by (Wheatland et al., 2000) which has been used to compute 
force-free fields. 



3. Convergence Tests 

Since analytic truly 3D MHS equilibria are not available, we use an 
analytic 2D MHS-equilibrium to test the newly developed code. The 
analytic equilibria are not meant to be a good representation of the 
solar corona and the tests are only carried out to check the convergence 
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of the newly developed code. We represent the magnetic field with help 
of the flux-function A(x, z) as 

B = VA x e y + By e y (12) 

and the MHS equations reduce to a Grad-Shafranov equation 

8 ( Bl{A)\ 
-AA = — (p(A^) + ^^\. (13) 



3.1. Equilibrium MHS-1 

As a first test we consider equilibria without gravity (P = P(A)) and 
choose P(A) oc A 2 and B y (A) oc A. The corresponding Grad-Shafranov 
equation is linear in A and can be solved analytically by a separation 
ansatz. It is convenient to define a function 

11(A) = c 2 A 2 = P(A) + S^-, (14) 
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P{A) = a c 2 A 2 , (15) 
By (A) = - oo) c A. (16) 

Configurations with c = correspond to potential fields, finite c and 
ao = to force-free equilibria, ao = 1 to equilibria with pressure gra- 
dient but without magnetic shear and finally < ao < 1 is the general 
case including both forces and magnetic shear. With this approach we 
get the solution of (13): 

oo 

A(x, z) = exp(— uttz/L) [a,k cos(kirx / L) + bk sin(fc7rx/L)] , (17) 
fc=i 

where v = \Jk 2 — c 2 , for c 2 < k 2 . The solutions of the Grad-Shafranov 
equation are invariant in one spatial coordinate (^ = 0). To test our 
3D-optimization code it would be more convenient to have equilibria 
varying in all three spatial directions. We construct such equilibria by 
rotating the solution of Grad Shafranov equation by an angle <\>\ around 
the z-axis and by (f>2 around the y-axis. As a result the solution varies 
in all of our three coordinate directions. The final equilibrium has the 
following free parameters: c, ao, a&, </>i, 02- As an example we choose 
(c = 0.8, ao = 0.5, ai = 1.0, 0,3 = — 0.8, au = bk = for all other k, 
^ = -0.057r,<fo = 0.15tt.) With help of the flux function (17), the 
equations for the pressure (15), shear field (16) and the magnetic field 
definition (12) we get the magnetic field B and the plasma pressure 
P. For an isothermal plasma we derive the density as p = -j^. We 
normalize the maximum normal magnetic field at z = to 300 Gauss 
= 0.03 T. 

As a first test, we want to reconstruct this equilibrium with our 
code. The code needs any 3D-vector field as start configuration for the 
iteration procedure and it is convenient to choose a potential magnetic 
field with respect to the photospheric line of sight magnetic field. The 
potential field can be easily computed for an observed magnetogram. 
For our model data a potential field is computed with the same pa- 
rameter set but c = 0. The boundary values for the iteration are in 
practical cases only known on the bottom plane. On the others plans 
they have to be iterated too using (11). In this test we fix the magnetic 
field however everywhere on the boundary to the value of the analytic 
solution to simplify the problem. For the force-free case we treated the 
side and top boundaries as unknowns in (Wiegelmann and Neukirch, 
2003) and showed how they can be iteratively determined by (11). 
This latter way of treating the boundary values makes the finding of a 
solution much more difficult. 

In figure 1 we show three-dimensional plots of selected field lines 
for this MHS-equilibrium. The colour coding of the bottom boundary 
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indicates the distribution of B n on that boundary. To test our code we 
try to reconstruct that equilibrium in the following way: 

— Inside the computational box we choose a potential magnetic field 
as start equilibrium. 

— We prescribe the plasma pressure in the box with the analytic 
solution. 

— We prescribe the vector magnetic field on the boundaries of the 
computational box. 

— We iterate for the magnetic field inside the computational box with 
(9). 

During the computations we calculate the quantities L, the absolute 
value of the force balance | J x B — VP\ (averaged over the numerical 
grid), the value of |V • B| (averaged over the numerical grid), and the 
difference between the numerical magnetic field and the known analyti- 
cal solution |B(i)— B ana | 2 /|B ana | 2 (averaged over the numerical grid) at 
each time step. In figure 2 we show the development of these quantities 
during the iteration process with logarithmic scaling. All quantities de- 
crease over several orders of magnitude during the optimization process 
and reach the level of the discretisation error. The discretisation error 
corresponds to the value of L , |J x B — VP — p V^| and |V • B| for 
the analytic solution computed on a numerical grid. In the upper half 
of table I we summarize the main result of the iteration process. The 
first row corresponds to the analytic solution computed on a grid and 
defines the discretisation error. The second row contains the values of 
L, the force balance and the relative error for the start configuration, 
where the interior points have been replaced by a potential field. The 
relatively large values of the three quantities are due to the deviation 
from the equilibrium. The next rows show how the three diagnostic 
quantities evolve during the iteration. After 10000 iteration-steps the 
discretisation error is reached for all values and the original equilibrium 
MHS-1 has been reconstructed. 

We use a Landweber iteration (see e.g. (Louis, 1989)) with some 
automatic control of the stepsize. The continuous form of equation (9) 
ensures a monotonously decreasing L. A monotonously decreasing L 
in the discretized form is ensured if the iteration step dt is sufficiently 
small. The code checks if L(t + dt) < L(t) after each time step and 
if the condition is not fulfilled, the iteration step is refused and dt 
is reduced by a factor of 2. After each successful iteration step we 
increase dt slowly by a factor of 1.01 to allow the time step to become 
as large as possible with respect to the stability condition. In principle 
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Table I. Details of runs to reconstruct MHS equilibria. 



Tlx X ^ ^2 


Step 


L 


Force - 


balance 


Relative Error 


7 — 9 7 


[nN 


m 6 ] 


40 x 40 x 20 


MHS-1 


0.0028 




0.19 


Reference 


Start 





140613 




1205 


0.27 




500 


74 




25.5 


0.025 




5000 


0.021 




0.4 


1.4 10~ 5 




10000 


0.0028 




0.19 


< 10~ 6 


40 x 40 x 200 


MHS-2 


4.3 




10.7 


Reference 


Start 





3.7 10 7 




15290 


0.14 




500 


101423 




2399 


0.017 




1000 


3393 




430 


9 10~ 4 




5000 


4.3 




10.8 


< 10" 6 



there are more sophisticated methods available to calculate an effective 
dt for each iteration step (see e.g. (Geiger and Kanzow, 1999)) but 
these methods have a huge numerical overhead and further numerical 
experiments will have to show which of these are favourable for our 
problem. 

3.2. Equilibrium MHS-2, Helmet Streamer 

As a second example we consider an equilibrium with gravity which has 
been used to model coronal helmet streamers (Wiegelmann et al., 1998). 
The method is based on the asymptotic expansion method (Schindler, 
1972) and corresponds to a nonlinear Grad-Shafranov equation. Here 
we choose for the terms in (13) 

P(A*) = a ex p("^) exp(cA), (18) 

By (A) = 2 exp Q A^j . (19) 

For simplicity we use a constant gravity f = 3 z, g = 270^ and a 
constant coronal temperature T = 3 10 6 K . Consequently we get the 
mass density from the plasma pressure as p = P/RT. The parameters 
correspond to a coronal pressure scale height 

_ k B T _ RT 

Ho ~ hv^ " — (20) 
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Log ( L) Log (UxB — VP I / [ n N/m-3 ] ) 




O 2000 4000 6000 SOOO O 2000 4000 6000 SOOO 

Steps Steps 



Figure 2. MHS-1: Evolution of L, force balance |j x B - VP|, |V ■ B| and the 
difference between the numerical magnetic field and the known analytical solution 
|B(t) — B ana | 2 /|B ana | 2 . All quantities are averaged over the numerical grid. 



of 93 Mm f» 0.13 solar radii. With help of the method of asymptotic 
expansion we find an analytic solution of (13): 




We choose c = 0.05 Mm 1 and A = 0.001 Mm 1 for the free parameters 
and compute the solution on a grid of nx = ny = 40, nz = 200. With 



twbi03.tex; 6/02/2008; 17:41; p. 10 



Magnetic modelling and tomography 



11 




Figure 3. MHS 2: A projection of some field lines for MHS-2. The back ground 
colours correspond to the logarithm of the electron number density p^j. Please 
note the different scale in x and z. 

help of the flux function (21), the equations for the pressure (18), shear 
field (19) and the magnetic field definition (12) we get the magnetic 
field B and the plasma pressure P. For an isothermal plasma we derive 
the number density as N = k ^ T - Let us remark that the quantity N 
is what we will get for real data from coronal measurements after the 
tomographic reconstruction. In principle a non-constant temperature 
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Log(L) Log ( I J x B — VP — p V* I ) 




Steps Steps 

Figure 4. MHS-2: Evolution of L, force balance |j x B - VP - pV*|, | V ■ B| and the 
difference between the numerical magnetic field and the known analytical solution 
B(t) — B ana | 2 /|B ana | 2 . All quantities are averaged over the numerical grid. 



T can also be used, e.g. from a standard atmosphere model. Figure 
3 shows a projection of magnetic field lines and the electron density 
structure as background for the helmet streamer configuration MHS-2. 

The solution is invariant in y and we rotate the solution around 
the z-axis with an angle 4> = n/10 which is useful for testing our 
reconstruction code (All derivatives appear) . As start magnetic field we 
choose a Harris-sheet, where the magnetic field has only one component 
B z = -V2 tanh(-^ x) which is invariant in z and y. We apply our 
code for the reconstruction of this helmet streamer equilibrium MHS-2 
in the following way: 
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— Inside the computational box we choose a Harris-sheet magnetic 
field as start equilibrium. 

— We prescribe the electron number density N in the box with the 
analytic solution. 

— Under the assumption of a constant coronal temperature and a 
constant gravity we calculate the plasma pressure, density and 
compute u = — /xo(VP + pV*) on the grid. 

— We prescribe the vector magnetic field on the boundaries of the 
computational box. Similar as in the previous example we use the 
analytic solution to fix the magnetic field on all boundaries for 
simplicity. 

— We iterate for the magnetic field inside the computational box with 
(9). 

Similar as for MHS-1 we diagnose the quantities L, the absolute value 
of the force balance | J x B — VP — p V'I'I (averaged over the numerical 
grid), the value of |V • B| (averaged over the numerical grid), and the 
difference between the numerical magnetic field and the known analyti- 
cal solution |B(i)— B ana | 2 /|B ana | 2 (averaged over the numerical grid) at 
each time step. In figure 4 we show the development of these quantities 
during the iteration process with logarithmic scaling. All quantities de- 
crease over several orders of magnitude during the optimization process 
and reach the level of the discretisation error. The discretisation error 
corresponds to the value of L , | J x B — VP — p V^| and | V • B| for the 
analytic solution computed on a numerical grid. Please note that the 
discretisation error for MHS-2 is significantly larger as for MHS-1 due 
to the nature of nonlinear analytic solution. In the lower half of table 
I we summarize the main result of the iteration process. 

The first row corresponds to the analytic solution computed on a grid 
and defines the discretisation error. The second row contains the values 
of L, the force balance and the relative error for the start configuration, 
where the interior points have been replaced by a Harris-sheet magnetic 
field. The large values of the three quantities are due to the deviation 
from the equilibrium. The next rows show how the three diagnostic 
quantities evolve during the iteration. After 5000 iteration-steps the 
discretisation error is reached for all values and the original helmet 
streamer-configuration MHS-2 has been reconstructed. 
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3.3. FORCE-FREE AND MHS RECONSTRUCTION FOR DIFFERENT (3 

It is generally assumed that the magnetic pressure in the lower corona 
is much larger than the plasma pressure leading to (3 <C 1. For longer 
structures like helmet streamers the plasma (3 increases significantly. 
It is generally assumed that the effects of plasma pressure and gravity 
can be neglected for low (3 plasmas leading to a nearly force-free state. 
Let us remark that one can construct high (3 force-free equilibria by 
adding a homogeneous plasma pressure VP = (or a barometric 
density distribution VP = —p V^ for configurations with gravity) to 
any exact force-free configuration. Here we do not study such singular 
cases, but more realistic configurations where VP oc j, where I is a 
typical length scale of the problem. The equilibrium MHS-1 allows us 
to compute configurations with different plasma (3 by prescribing the 
parameter ao, where ao = corresponds to an exact force-free state. 
We use our code to investigate how well a magnetic field configuration 
can be reconstructed by the force-free approach. We start all recon- 
struction runs with a potential field solution c = 0.0 and the magnetic 
field boundary conditions extracted from the exact solution (artificial 
vector magnetograms). The configurations in the left hand side in table 
II have been reconstructed by the assumption of a force-free configu- 
ration and the configurations in the right hand side in table II are 
reconstructed as MHS-equilibria. The force-free reconstruction needs 
the boundary magnetic field data as input and the MHS-reconstruction 
additional requires the plasma-density structure, which we extract here 
from the analytical solution. This corresponds to artificial tomographic 
information. 

We diagnose the quantities L (Lff for ForceFree and Lmhs f° r 
MagnetoHydroStatic) and the deviation from the the analytic solution 
(ErrorFF for ForceFree and ErrorMHS f° r MagnetoHydroStatic) sim- 
ilarly as described in the previous sections. For MHS-reconstruction 
our code finds the magnetic field structure for all configurations with 
an error corresponding to the discretisation error. If we restrict our 
code to an exact force-free reconstruction we still get a considerable 
good agreement with the exact solution for a plasma (3 of less than 
10~ 3 . For higher values of (3 both the value of L (where L = cor- 
responds to an exact force-free state) and the error in the magnetic 
field increases. Consequently our code finds the expected result that 
the effect of plasma pressure is neglible for low (3 configurations. The 
result also shows that a direct consideration of tomographic information 
regarding the electron density is only useful for finite (3 plasmas. For 
low (3 plasmas the magnetic field structure is practically not influenced 
by the plasma density distribution. Let us remark that it is still possible 
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to consider some indirect information provided by the plasma density 
for low (3 plasmas, e.g. the fact that the density gradient parallel to the 
magnetic field is much lower than the density gradient perpendicular to 
the magnetic field. Consequently the magnetic field lines are outlined by 
the emitting plasma. This allows to consider stereoscopic information 
for the reconstruction of low (3 plasmas (Wiegelmann and Neukirch, 
2002). 

3.4. Speed of the method 

The speed of our code is approximately proportional to N 5 (N is the 
number of points for one side of the computational box) similar as found 
by (Wheatland et al., 2000) for the force- free case. This iV 5 dependence 
looks discouraging for the reconstruction of large boxes. We undertake 
some rough estimations if the method is practical for modern vector 
magnetographs. For a grid of N=40 a reconstruction takes about 5 min 
on a 4 processor computer. The IVM vector magnetograph in Hawaii 
has a resolution of N=256 pixel. Consequently a reconstruction with full 
IVM-resolution would take approximately 5min • (256/40) 5 « 35days. 
If only the half IVM-resolution is used the reconstruction time would 
be 5min • (128/40) 5 ~ 28hours, which seems to be acceptable. We are 
optimistic that an improved numerical scheme (e.g. using conjugated 
gradients or multi-grid methods), a massive parallelization (Using 16- 
32 processors instead of 4. 1 ) and increasing computer speed will speed 
up the reconstruction time significantly. The computing time for the 
optimization approach seems to be high, but comparisons of the opti- 
mization method with classical MHD relaxation (for the force- free case) 
have shown that the optimization method is more effective (Wiegel- 
mann and Neukirch, 2003). Direct extrapolation methods (e.g. (Wu, 
Chang, and Hagyard, 1985)) are much faster than iterative methods 
but known to become unstable with increasing coronal height. 

4. Using coronal magnetic field information as a constraint 

for tomography 

Until now, we used information regarding the coronal density and 
pressure structure as given. In principal we could consider density and 
pressure in functional (4) as additional variables to be optimized just 
as B. In this case, however, the problem of minimizing L would be 



1 The method seems to parallelize quite well. On 4 processors the reconstruction 
is about 3 times faster than on one processor. 



twbi03.tex; 6/02/2008; 17:41; p. 15 



16 



Wicgclmann and Inhcster 



Table II. force-free and MHS-reconstruction. The first column 
contains the plasma /?, the second column oo, the third column 
the final value of L for a force-free reconstruction, the fourth 
column the error in the magnetic field structure compared to 
the analytic solution, the fifth column the final L-value for a 
MHS-reconstruction and finally the sixth column the error in 
the magnetic field structure compared to the analytic solution 
for the MHS-reconstruction. All runs have been computed with 
MHS-1 on a grid n x = n y = 40, n z — 20 for 5000 iteration 
steps. 



Plasma (5 


fflo 


Lfp 


ErrorFp 


Lmhs 


Error M HS 








0.027 


1.5 10~ 5 






10" 4 


0.00045 


0.028 


1.5 10~ 5 


0.028 


1.5 10~ 5 


10~ 3 


0.0045 


0.04 


1.5 10~ 4 


0.028 


1.5 10~ 5 


10~ 2 


0.045 


1.26 


4.9 10~ 4 


0.028 


1.5 10~ 5 


HT 1 


0.41 


112 


6.0 10 -3 


0.027 


1.7 10~ 5 


0.2 


0.74 


807 


0.05 


0.015 


1.0 10~ 5 


0.3 


1.0 


1443 


0.2 


0.007 


2.0 10~ 6 



hopelessly underdetermined even if the boundary values of B, p and P 
were given at the Sun's surface. 

If the magnetic field was known on the other hand, we could use (4) 
and assume a temperature variation to determine a consistent density 
and pressure. Even though we will find immediately that this approach 
is doomed to fail we here mention for completeness the optimization 
equations which can be derived from (4) 

p = mN, P = k B TN 

and we get (see the appendix B) 

1 dL f TT 8N „ f T dN l2 . A X 

H = p m n a ■ - no k B T V • ft a , (25) 
I = HO k B T ft a ■ n. (26) 

(See (8) for the definition of fi a and fib-) 2 Physically speaking, the 
solution to these equations yield a pressure which exactly balances ex- 



2 Let us remark that the general form of ^ will vary both the density distribution 
and the magnetic field 

1 dL f dB „ rr dN , 3 f dB „ r dN , 2 

2lE = ~J v lH- F + H ln dx ~J s lH- G + I lH dx - 
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cessive magnetic field forces onto the plasma. However, due to the small 
value of j3 in the corona, already a small relative error in the assumed 
magnetic field will result in residual forces which need a huge plasma 
pressure to be balanced so that a small relative error in B will produce 
a large relative error in N. We obviously need another approach to get 
hold of a decent estimate of the coronal density distribution and this 
must be is based on additional observations. 

The coronal electron density can be observed by coronagraphs. Un- 
fortunately, coronagraphs yield only integrated column densities along 
the line-of-sight and the 3D density distribution itself has to be recon- 
structed from these measurements by means of a tomographic inversion 
(Zidowitz, 1999; Frazin, 2000; Frazin, 2002). This inversion process has 
besides its intrinsic ill-posednes additional problems to cope with if 
applied to coronagraph data: 

• non-stationarity of the coronal density structures, 

• incomplete data due to the occultation of the image centers (exterior 
tomography problem), 

• non-ideal viewing geometry due to the tilt of the Sun's axis with 
respect to the ecliptic. 

As a result of all these problems, the spatial resolution which ultimately 
can be achieved with the reconstruction is limited. The conventional 
procedure to obtain a reliable solution is to minimize an expression of 
the following form: 

G(N) =£l<f- l P ,m\ 2 +» J \R(N)\ 2 d 3 x, (27) 
P>* v 

where I^f is the observed intensity in pixel p of image i and X Pi j is the 
respective simulated intensity which is calculated from a density N as 
a line-of-sight integral 

I p ,i = J Ndt (28) 

Here, C Pi j is the beam from pixel p of image i, i.e. the location of points £ 
V which project onto the respective pixel. We omit here modifications of 
(28) due to the scattering geometry and the scattering crossection of the 
electrons which lead to slight variations of the integrand in (28). In (27), 
R is a regularization function to be specified below. The primary aim is 
to find a density model N which makes the first term vanish. In this case 
N is compatible with the observations. However, due to measurements 
errors, inconsistencies of the observations mentioned above and its pos- 
sible insensitivity to certain density structures, it does not make sense 
to minimize the first term alone below a level approximately given by 
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the measurement error variance. To stabilize the model reconstruction 
on density structures to which the observations are insensitive, the 
regularization term is added with a regularization parameter fi tuned 
so that the first term is approximately brought down to the observation 
error variance when the complete expression G is minimal. 

Since I(N) is basically an integration operator, R(N) very often is 
chosen as a differential operator, e.g., 

R(N) = V 2 N (29) 

As a result, while the integration makes I{N) insensitive to small scale 
structures, R(N) responds to these but has little effect on the large 
scales. Therefore, when minimizing G, the large scales are shaped by 
the first term to comply with the observations, while the small scales are 
kept smooth by the regularization term in G. The transition between 
large and small scales is mainly determined by the weight fi of the 
regularization in G which must be enhanced the worse the quality of the 
observations. In this sense, the spatial resolution of the model density 
TV which can be achieved depends largely on errors, inconsistencies and 
gaps of the data. 

But the power of the regularization term goes beyond suppressing 
unwanted small scale noise in the reconstructed model N. Any ad- 
ditional physical constraint for N can be included here just like the 
V • B = constraint was added to the force balance condition in (4). 
One obvious constraint for the density is its positivity, which can be 
taken account of in R by so-called barrier functions (see (Prazin, 2000)). 
For a more stabilizing constraint we may return to (24). It was derived 
so as to minimize the square of f2 a which is proportional to the local 
force balance. Note that fib has no dependence on iV and therefore 
is not varied here. The argument which led us to discard (24) as a 
starting point for an iteration for the density mainly applies to the force 
components tt a ± oc (V x B) x B — /io(Vj_P + pVj_^) perpendicular 
the magnetic field which are dominated by the magnetic term to lowest 
order. We therefore choose a regularization term which at least takes 
care of the field-aligned force balance in fi 0j || oc — /xo(V||P + pV^). 
This is achieved by 

R(N) = -L(B • V)P + p(B ■ V)* = (B • V)N - ^N, (30) 
k B T H 

where again an isothermal plasma is assumed and Hq is the pressure 
scale height of the solar corona as in (20). An extension to a varying 
temperature is straight forward if it is given, however we cannot solve 
for an unknown T unless we make use of additional observations. 
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Figure 5. A reconstruction of completely noisy images with the regularization term 
(30). The visible structures are imprints of the dipole magnetic field assumed in the 
regularization term. 

The regularization term (30) effectively smoothes the density out 
along the magnetic field lines and thereby takes account of the fact that 
the transport coefficients in a magnetized plasma are much larger along 
the magnetic field than in perpendicular direction. In perpendicular 
direction to B the density may have gradients as sharp as our model 
resolution allows without changing the value of G. 

The magnetic field assumed in the test calculations below was a sim- 
ple dipole field. The effect of the above regularization operator becomes 
particularly visible if instead of meaningful data, we assume that I obs is 
pure noise. The magnetic field then is the only real information in the 
inversion process and it becomes directly visible in the reconstruction 
results (see Fig 5). 

In appendix C we derive an expression for an iterative descent step 
analogous to (24) but which is preferable to (24) because it includes 
the additional observations to stabilize the reconstruction 




(31) 



v 



s 



where 



2^6 Cp , i (x){l p , i (N)-lZf) 



2a4(B- V)R(N) + ^R(N)], 



(32) 
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J(x) = 2n(B-n)R(N) (33) 

and 5c p { is 1 inside the beam emanating from pixel p of image i and 
otherwise (see appendix C). 

In order to test the feasibility of our scheme, we minimize (31) for a 
2D test model by means of a conjugate gradient iteration. In table III 
we compare the action of this operator with more conventional means 
of regularization for a two-dimensional reconstruction. Here, the model 
and data errors of a reconstructed model density N are defined as 

x Pll (A0l 2 , 

N\ 2 d 3 x. 

Here, iV ana is the analytic density model used to obtain the data I obs . 
The reconstruction a) was obtained without explicit regularization (i.e. 
with fi = 0). Instead, the iteration was stopped after 13 iterations when 
the model error was minimal. Subsequent iterations further decrease 
the data error but enhance the model error (which for real data is 
not known), a phenomenon which is known as semiconvergence. For 
reconstructions b) and d) the conventional regularization operator (29), 
for c) and e) the operator (30) was used, however with the scale height 
Ho set to oo. 

In cases d) and e), the regularization parameter fi and the number 
of iterations where optimized to achieve the best agreement with the 
original model. We show these results only to demonstrate how close a 
reconstruction can come to the true solution in principle. In practical 
cases, however, the true density distribution is not known and the 
optimum solution has to be sought based only on the values for the 
data error and the regularization term. In cases b) and c) these values 
have been used in an L-curve search for the optimum solution (Hansen 
and O'Leary, 1993). 

In Fig. 6 we show the models associated with the test inversions in 
table III. The upper left shows the original model simulating a coronal 
loop on the western limb and a streamer on the eastern limb. This 
model was used to calculate the artificial data used as input for the 
reconstruction after some noise was added. The standard deviation of 
the noise was 3-10 -2 times the square root of the local data intensity 
when the maximum data intensity is normalized to unity. 

In terms of the model error, (30) yields slightly better results than 
(29). The major improvement comes from the region close to the oc- 
culter. In this region conventional tomography can only yield a limited 



data error = -^T^ 8 

p,i 

model error = - J \N ana 
v 
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Table III. Comparison of different regularization methods. The model and 
data errors are defined in the text. The results a) to c) correspond to the 
cases in Fig. 6. Case a) uses no regularization at all and has a minimum 
model error after 14 iteration steps. Cases b) and d) use the isotropic 
regularzation as in (29), cases c) and e) the isotropic regularzation as in 
(29). Iteration steps and n are optimized in b) and c) from an L-curve 
search, in cases d) and e) for a minimum model error.) 



Regularization method 


optimal iteration 
steps 


model error 


data error 


a) stopping rule, R=0 


14 


12.4 


6.46 10~ 5 


b) R=(29), ^=0.01082 

c) i?=(30), ^=0.00787 


102 
103 


1.65 
1.21 


13.9 10~ 5 
12.5 10~ 5 


d) i?=(29), ^=0.0050 

e) R=(30), ^=0.0065 


53 
120 


1.21 
1.13 


12.2 10~ 5 
12.0 10" 5 



resolution because close to the occulter the structures are only seen in 
few observations. On the other hand, this is the area where we observe 
the strongest gradients in the density structures and where spatial 
resolution is needed most. This fundamental lack of resolution can only 
be overcome if new information (either observations or assumptions) is 
fed into the inversion process. The new regularization operator (30) 
contains this additional information in form of the local magnetic field 
direction. 

The price we have to pay is a more lengthy computation as the num- 
ber of iterations increases (see Table III). Another problem which might 
occur are spurious field-aligned density structures in the reconstruction 
which add up to zero in the tomographic projections. Formally, (30) 
has a nullspace which asymptotically comes close to the nullspace of 
the 2. In practical cases, however, the discretization error in the dif- 
ferentiation is sufficient to give (30) some isotropic component so that 
even exactly field-aligned structures yield a small non-zero contribution 
in a discretized operator (30). 



5. Conclusions 

In this paper we undertook a first step towards the inclusion of to- 
mographic information into the reconstruction of coronal magnetic 
fields and a first step towards the inclusion of coronal magnetic field 
information into the tomographic inversion procedure. Until now we 
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Figure 6. Comparison of different regularization methods for the reconstruction of 
the original model -/V ana in the upper left. The results a) to c) correspond to the 
cases in table III. 



considered the reconstruction of the magnetic field from its boundary 
values with an optimization code when the density structure is given 
and the tomographic reconstruction of the coronal density distribution 
from coronagraph data under the constraint of a given magnetic field. 
As neither N and nor B are known a priori in the solar corona, we 
rather have to find a way to consistently reconstruct both quantities 
from the observations simultaneously without the assumption that one 
of the quantities is given. 

In fact we observe that the two approaches discussed in the previous 
sections are not only formally closely related but can be derived from 
a single variational problem if the expression for L is slightly modified 
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and the factor B 2 



in the integrand is omitted 



L(B, N) = f |(V x B) x B - /i (VP + /)Vf )| 2 + B 2 |V • B| 2 ] d 3 



Jv 



X 



(34) 



p,i 



A variation with respect to B obviously leads to an iteration scheme 
similar to (5), except for the effect of the B~ 2 term in the integrand. 
We have tested the resulting scheme and found that it also converges 
towards the analytic solution from which the boundary conditions were 
taken, but the convergence speed was much slower than with (5). If L is 
varied with respect to N alone, we can ignore the terms which depend 
only on B and we obtain (31) again if we discard the perpendicular force 
balance in fi a . This suggests that the individual reconstruction prob- 
lems for B and N are just two projections of a unique reconstruction 
problem. In this case we could apply both algorithms simultaneously 
and replace N in the algorithm for B and B in the algorithm for N by 
the respective iterate and the problem as a whole should converge as 
they do individually. A test of this hypothesis will be attempted in the 
future. The code is planed for use within the STEREO mission. 
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Appendix 



A. Derivation of F and G in (5). 



(4) can also be written as 




(A.l) 



ft a = B~ 2 [(V x B) x B + u] 
ft b = B- 2 [(V-B) B] 
u = -/i (VP + pVtf) 



(A.2) 
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We vary L with respect to an iteration parameter t and get 

1 dL f „ d 

2 7ft 



= ^fi a - ^[(V x B) x B + u] d 3 x 
+ / v nb-|[(V-B)B]c?x 



(A.3) 



Our aim is now to use vector identities and Gauss law in such way 
that all terms contain a product with This will allow us to provide 
explicit evolution equations for B to minimize L. The third term has 
the correct form already. We expand the first and second term 



2 dt ~ Jv 8 



V x 



<9B\ 



+ 
+ 
+ 



ft, 



v 



dt J 
(V x B) x 



I x B 
dB 



dt 



d 3 x 
d 3 x 



Jv 



V 



b • 



d s x 



d 3 x 



v (n 2 a + nl)B.^d 3 x. 



(A.4) 



The fourth and fifth term have the correct form. We apply the vector 
identities a • (b x c) = b • (c x a) = c • (a x b) to the first and second 



term 



I dL f / <9B\ ,„ n . ,o 



+ 
+ 



^ — -(n a x (VxB)) d 3 x 
dB 

n b -B) v-^-d 3 * 

dB 



+ J> b (V.B)].-A 



(A.5) 



Term two, four and five have the correct form. We apply (V x a) • b = 
a • (V x b) + V • (a x b) to term 1 and ipV ■ a = -a • + V • (aifj) to 
term 3 

1 ^ = -/ y ^-[Vx(ftaXB)]d 3 X 



2 dt 
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+ 



dB 



(f2 a x B) x 



u • (O a x (V x B) d 3 x 
v dt 

-X v(nb - B) -f A 

(n b .B)— ci 3 x 
^ [n b (v • b)] • — d 3 x 



+ v- 



(A.6) 



Term one,three, four, six and seven have the correct form. We apply 
Gauss law to term two and five 

1 dL f dB 



2 dt 



-/A 



dt 



• [V x (O a x B)] d 6 x 



(fi a x B) x 



dB 

~dt 



d 2 x 



+ J ^ • ("a X (V X B)) d*X 

+ / 5 n(« b .B).f ^ 

+ |[ fib (V.B)].-^ 

- / (^ + ^)B 



,2 , o2^ dB d 3 x . 



Iv - at— (A - 7) 

Now all terms but the second have the correct form. We apply a • (b x 
c) = c • (a x b) to the second term 

1 dL f dB r _ , n ,„ 

=*2 * = -/^•[VxtSi.xB)]^ 

-/[fix [(n, x b)] ~ A 



/J 

— • (O a X (V X B)) d 3 X 

-I 

+ / sfl(SV B).f t 



dB 

v(n h ■ B) ■ — d 3 x 

V OT 

<9B 
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+ | [fM v.B)].f A 

/<9B 
(^ + ^)B- — d 3 x. (A.8) 

Now all terms have the correct form and we write them more compact 
as 

1 dL f dB „ o f dB „ l2 /a \ 

=>- — = -/ — - • F d 3 i - / — - • G d 2 x, A.9 

with 

F = V x (fi a x B) - fl a x (V x B) 

+v(n b • b) - n b (v • b) + (nj + ng) b, (a.io) 

G = n x (n a x B) - n{ft h ■ B). (A.ll) 



B. Derivation of H and I in (24). 

We vary L with respect to an iteration parameter t where the magnetic 
field is independent from t here. 

= f n a ~[(V xB)xB- M (VP + /9V*)]d 3 x (B.l) 
2 at at 

We write the pressure P and the mass density p as functions of the 
number density A 

- /lofeTVJV + mJVV^A (B.2) 

Our aim is to write all terms as products with ^ to derive evolution 
equations for N to minimize L. The first term vanishes here because B 
does not depend n t. We expand the second term 

■*5f--"X^j(£)'* 

-/x / fia-m — Wd'V (B.3) 
at 

The second term has the correct form. We apply a-X7ip = V-(a?/>)— ?/>V-a 
to the first term 

1 dL f „ / , m ajv\ o 
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dN 



+fi f V • ft a k B T — d 3 x 

J^n a -m — W d 3 x. (B.4) 



Term two and three have the correct form. We apply Gauss law to the 
first term 



1 dL f n dN j2 

2Tt = J s n-Cl a k B T — dx 

f dN , 

+Ho / V • ft a k B T — d 3 x 
Jv ot 

f dN o 

-HO / ft a • m — ci 3 x. (B.5) 
Jv ot 



Now all terms have the correct form and we write more compact 

1 dL f ON , 3 f T dN ~ 

=> = - / H — d 6 x - I — d 2 x 

= jtxom fi a • V* - /u fc B T V • fi a 
J = Mo /c B T fl a • n. (B.6) 



C. Derivation of H and I in (31). 

An essential advantage in the derivation of the variational derivative 
of G in (27) is the fact that G is convex quadratic expression. While 
jobs j g j us ^. a data vector, X is a linear operator from the model space 
{N(V)} into the data space {I p ,i}- 

Z P) i(N) = j N(x)d£ : model space — ► data space. (CI) 

The beam was defined in the main text as C Pt i = {x £ V | x projects 
onto pixel p for the view direction of image i }. The adjoint of (C.l) is 

i7(x, I) = ^2 $c p i ( x ) Ip,i '■ data space — ► model space, (C.2) 

p,i 

where 5c pi (x) = 1 for x £ C Pt i and Sc pi = else. The adjointness can 
easily be checked by insertion of the respective definitions into 

Y,I P ,a P A N ) = J J(*,I)N(x)d 3 x (C.3) 
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Now minimizing (27) 

G(N) = £ - 1 P ,(N)\ 2 +nj R(Nf S 

with (30) 



x 



V 



R(N) = _L(B • V)P + p(B ■ V)* = (B • V)N - 



yields 



p,i 

+ 2fi J R 

v 



dt 

, A „.dN B-rdN 

(BV) iir-ifri« 



d 3 x 



2jj(x,l P:i (N)-I°f] 
v 



dN 
~dt 



2/z 



v 



(B-V)R + R 



B r 



Ho 



~dt dx 



+ 2n J R(B-n)^d 2 x 

s 



V s 



where n is the unit normal on the surface S and 

#(x) = 2^s Cp ^)(i P ,m-i P t) 



p,l 



2p [(B-V)R(N) + B -^R(N)], 
-Ho 



/(x) = 2/i (B • h)R(N). 



(C.4) 



(C.5) 



(C.6) 
(C.7) 
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